TIME DEPENDENCE OF THE MAGNETIC 
FIELD IN A RECTANGULAR TOROID 


by 
44S 
GLEN LEROY SHURTZ 


B.S.M.E., Kansas State University, 1963 
B.S.E.E., Kansas State University, 196 


A MASTER'S THESIS 
submitted in partial fulfillment of the 


requirements for the degree 
MASTER OF SCIENCE 


Department of Electrical Engineering 


KANSAS STATE UNIVERSITY 
Manhattan, Kansas 


1966 


Approved by: 


BL Meh ety wk 
SALE EIEE CTE 


TABLE OF CONTENTS 


‘§ $¢2 


1.0 
2.0 


: AcUmM#t- 7 


EEEPODUCTION. .. o: ble 0) &y.0! 0b a 6 6,0 (es Nelebl et es 


DERIVATION OF THE HYSTERETIC DIFFUSION EQUATION 
FOR A MAGNETIC BOUNDARY VALUE PROBLEM ..... 


2oa. IntroductiGis tc + + Geral els) ol) ee >) 6.6 
2.2 Magnetic Circuit Concepts .....e++-eee-s. 
2<3.. Eddy Curregmel is « «ie » svindelle oe 14.648 "18 
2.4 Flux Distribution in the Rectangular Toroid . 


2.5 Derivation of the Hysteretic Diffusion 
Equa tion , . . . . * . * . * . . * * . . * 


Maxwell's Equations ...+«-s+ssb#e+e#e 6» 
2.6 Approximation of the B-H Curve for U(H) ... 
2.7 Hysteretic Diffusion Equation ..... =. « « 
NUMERICAL SOLUTIONS WITH THE MODIFIED EULER METHOD 


3.1 Euler Method for Solving an Ordinary 
Differential Equation . ....+#-e« + s8 « » 


3.2 Modified Euler Method for Solving an Ordinary 
Differential Equation « . . « « « « « 6 6 « 


3.3 Numerical Solution of the Hysteretic Diffu- 
sion Equation with the Modified Euler 
Method . ° . . . ’ . . * . . . ’ * . . . * 

3.4 Conditions for Convergence .....».++s « 


3.5 Relationship Between At and Other Parameters 
to Insure Convergence . .... ++ « « « s 


PROGRAM FOR SOLUTION OF THE HYSTERETIC DIFFUSION 
EQUAT ION . . . *. . * . . J * . . * . * . . . . * 


Uodes LACPOGUNELOM « « « te af wale sen enw’ Guat ve. el ereanee 
if, 2: Main ‘BrogremiEMBX §. Aa ch a to: Shicus us SOnw 6 ee 
h. 2A . Comments (MG90000-MG90106) . . . . . . . 


4.2B. Input Data (MG90122-MG90138) and 
(MG90262) . . . . . . . . . . . . . . 


ii 


iii 


4.2C. Evaluation of the Hysteresis Loop 
oxi naa aay Formula iid al 
MG90250) .. o eww 50 


4.2D. Calculation of Constants and Page 
Headitic Rowbane ws: +. ella a ee ruta ae oe 


l..2E. Initial Conditions Print-out Routine 


(MG9I0 72 SMGGGC TIE) iahu «ohn. dee ck 0) a: ae ey SL 

4.2F. Programming the Modified Euler Process 

for the Hysteretic Diffusion Equa- 

tion (MG90876-MG91060) ....... Sil 
4.2G. Output Data (MG91078-MG9110) . iak wire 56 
l..3. “Siieemogram NEUE gyuts Pile sl oetda seals ss i a] ee ee ee 
Leal: 9 Sige pram “FORGE! yates ban Sos. on ad aatibas Nac) at, rooney sna 
lio... Subpmor ream Mealmt ts, myc ke gee gy yc yolks Wy big t's PEL 
5.0 USING THE PROGRAM ON A SAMPLE PROBLEM ....... 6h 
| Sie tepaTine the. Tait, EGA sane meats cis) oo te eke PROIE 
Beers, OUT SDA Uae oo. ate Mo deone ee eal meets SS ce oe vance 
5 ahiveOubpub Tapes 6 (Papen: (l-oe) ass s le oe wom 
B.cks OubputsTapes>’ (Paces:.63692) o.ts +) ss «sa Oe 
5,20. Outeut Vapi picrares os=200)k Me) sos lo Meee 
5.3 Discussion of Output Data and Results .'.. . 102 
ACK NOWLEDGMEINGS .. . egmaeieases VSP MEE aia'ol sl aus) Raumtislewk sos 6 cemtLOu 
RERERENCH Oise cts. .s Ma? gy eile t ot a> art so cb Site? wl.d oe! ak) eS 
APPONDICES © fir /-c-\s) sia > Saal ine okie ke” «is ts eeale LOO 


1.0 INTRODUCTION 


The purpose of this thesis is to formulete an equation 
which describes the time dependence of the flux distribution in 
& rectangular toroid of ferromagnetic material subject to given 
boundary conditions and to provide a numerical solution to the 
equation. This is of interest because it allows one to predict 
the time necessary to release a no-work magnet in an electro- 
mechanical system or define the switching time of a toroid used 
in the core plane of many modern digital computers. An analyt- 
ical solution also provides a convenient tool for optimizing 
the many problem variables or studying the effect of changing 
one or more of the variables on the magnetic performance of a 
given electromagnetic circuit. 

Development of theory and assumptions necessary to derive 
such an equation for the special case considered, i.e., for the 
toroid of rectangular cross section as illustrated in Fig. 7, 
or the electromechanical system as illustrated in Fig. 2, are 
given in section 2.0. The resulting equation obtained by 
manipulation of Maxwell's equations was found to resemble the 
diffusion equation when subject to the assumptions required for 
the problem of interest. The Hysteretic Diffusion Equation, 
equation (2-28), expresses the desired relationship for specify- 
ing the magnetic field intensity H as a function of position 
in the core and time. It should be noted that H is written as 


simply H for convenience and represents H(x,y,z,t). 


OH (co +H)*[ d°H dn 
(2-28) a 


where H is the magnetic field intensity (z component) 
Oo is the material conductivity 
Cy is a constant specifying material properties 


C5 is a constant specifying material properties. 


A discussion of the effects of eddy currents and concepts 
used to confirm the existence of a flux distribution in the cross 
section during a transient condition are given to aid the reader 
in understanding the problem. 

The task of evaluating H was accomplished by using the 
Modified Euler method of numerical integration in which the rec- 
tangular cross section was thought of as being divided into a 
grid. Each grid represented a toroid with a time-variant field 
intensity. This allowed one to express the Hysteretic Diffusion 
Equation in its finite difference form and approximate OH/dt 
by numerical techniques. Once this quantity was known, Euler's 
method was used to predict an approximate value for the field 
intensity of each grid area at a time At later, thus an 
approximate solution was obtained for H as a function of time 
for each grid area of the cross section. Values of flux density 
B, were obtained by substitution of H in an equation approxi- 
mating the B-H relationship for the given ferromagnetic material. 

One must remember the numerical process yields only an 
approximation to the actual values. Accuracy increases as the 


grid size decreases; however, the process becomes very slow for 


very small grid sizes because the maximum time increment allow- 
able to insure convergence of the Modified Euler process is 


dependent on many of the problem variables as shown by equation 


(3-44). 


2 
O-C,Coh 


(3-44) 0 <A} 
2(Co Hina x) 


Consideration of conditions for convergence and accuracy 
of the process are given in section 3.4. 

The last two sections of the thesis describe the program 
and various input variables needed to execute the process. 
Each input variable and its function in the program is discussed 
in section 4.2B. A sample problem with a typical input data 
set and the corresponding output data is given in section 5.0. 
Results of the data indicate a flux distribution and flux decay 


as predicted by sections 2.3 and 2.l. 


2.0 DERIVATION OF THE HYSTERETIC DIFFUSION 
EQUATION FOR A MAGNETIC BOUNDARY 
VALUE PROBLEM 


2.l1 Introduction 


A problem of interest to many is the solution of the dif- 
fusion equation sometimes called the heat equation. It is use- 
ful in describing the temperature distribution in an iron bar 
as a function of time, or in the magnetic case, the time depend- 
ence of the flux distribution in a toroid of magnetic material 
subject to given boundary conditions. It is the latter case 
which will be considered in detail in this thesis. The solution 
of this particular boundary value problem is important in pre- 
dicting performance objectives and analysis of many components 
in modern digital computers as well as being of use for analyt- 
ical evaluation of heat transfer characteristics and other 
problems in many other fields described by the equation. For 
' example, the memory of most modern computers consists of core 
planes in which rectangular toroids are affixed at the junction 
of two write windings as illustrated in Fig. l. 

To magnetize the toroid one-half of the current needed must 
flow in the proper direction through both write wires threaded 
through the core; thus the core is magnetized in one particular 
direction which designates a "1" bit and the opposite direction 
which designates a "0" bit. During a read out of the memory, 
the direction of magnetization determines the direction of cur- 


rent flow induced in a read winding and in turn determines the 


1/2 Current Selected core 


Fig. 1. Section of a memory core plane. 
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Fig. 2. Electromechanical device (no-work magnet). 


presence of a "0" or "1" bit. One can see that it may be desir- 
able to know the time to accomplish flux reversal during the 
write time. Furthermore, an analytical method of determining 
the reversal time for various cores of different sizes, shapes, 
materials, etc., would be convenient since this time could be a 
Significant factor in determining the speed of the computer. 

The solution could also provide information useful for im- 
provements in the design of present-day input-output equipment 
for computation systems. Many high-speed card punches and 
printers use electromechanical devices in various punching and 
printing techniques, all of which require elaborate sequential 
timing and mechanical movement to accomplish the desired result. 

| Let us consider the electromechanical device in Fig. 2. 
Suppose a voltage E causes a current I to flow through the 
coil which produces a flux % sufficient to overcome the force 
kx and hold the core and armature together. When I is re- 
moved (i.e., E = 0), the holding force is removed and the spring 
force allows the hammer to impact the stop. An extension of 
this principle is used for modern printers which can print at 
the rate of twelve hundred lines per minute with one hundred 
twenty characters per line. The time required to release the 
magnet becomes of prime importance since sequential timing and 
logic circuits required to release the armature cause a release 
and hold operation cycle to occur at very high rates of speed. 
A typical release time might be one millisecond. This same 
mechanism is used for obtaining punched cards and the same dis- 


cussion could apply. 


Due to effects of mechanical inertia, eddy currents gener- 
ated by the rapidly changing boundary conditions on the magnetic 
circuit, and other factors, it may be desirable to provide an 
analytical solution which would account for changes in the 
parameters affecting the electromagnetic performance of an elec- 
tromechanical system as outlined in the previous discussion. 
Assuming that the junction of the armature and core assembly of 
Fig. 2 does not provide an additional reluctance to impede the 
flow of flux across it other than that of the material itself, 
one can use the results of the following procedure to gain in- 
sight on effects that parameter changes produce on the electro- 


magnetic performance. 
2.2 Magnetic Circuit Concepts 


Magnetic circuit considerations are closely anslagous to 
those of resistive electrical circuits; however, the cause snd 
effect relationship in the magnetic case is nonlinear, i.e., the 
reluctance of a d-c magnetic circuit depends on the flux in the 
circuit, while for the d-c electrical case resistance is rela- 
tively unaffected by the amount of current. 

If one considers a toroid of magnetic material with a coil 
of wire wound tightly and distributed uniformly around it, a 
magnetic circuit problem is encountered. (See Fig. 3.) 

_ The voltage E produces a flux g in the megnetic material. 
The flux lines are perpendicular to a cross section of the toroid 


and should be uniformly distributed over the cross section in 


Pi te Site 


Fig. 3. Toroidal magnetic circuit. 
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Hysteresis loops and magnetization curve. 


the steady-state case. To insure uniform boundary conditions 
produced by the current I, assume the cross-sectional diameter 
of the ring is very small when compared to the inside and out- 
side diameters of the ring. Under these conditions the follow- 


ing relationships hold: 


(2-1) $= NI = Hl = RG 
(2-2) B = pH = @/A 

(2-3) R = pl/A 

(2-)) Force (F) = B°A/2u9 


The above also assumes that no leakage occurs and that one 
is only considering the steady-state solution; however, the 
transient operation is of major importance when an electromag- 
net or rectangular toroid is used in the applications as out- 
lined in section 2.1. 

If one desires to know the value of the flux at several 
discrete points in time during a build-up or decay cycle, one 
must use a different method of analysis than presented pre- 
viously. During the transient state a relationship exists be- 
tween flux and magnetic field intensity as illustrated by the 
hysteresis curve of Fig. 4. This curve specifies the magnetic 
properties of the ferromagnetic material. 

One can also see that during a flux build-up or decay re- 
sidual magnetism B,, and coercive force H,, depend on the orig- 
inal values of the flux density B and magnetic field intensity 
H. A typical ferrite material in the memory of a computer ex- 


hibits what is known as a square-loop property. This means that 
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the material's B-H characteristic has a square or rectangular 
hysteresis loop, which allows an immediate change in the direc- 
tion of magnetization once the proper value of H is present. 
In all cases the value of B for a corresponding H follows 
the upper curve if in the decay portion of the cycle and follows 
the lower curves if in the build-up portion of the cycle. 
Furthermore, during either build-up or decay, a flux distribu- 
tion exists across the cross section of the magnetic core due to 
eddy currents generated from the changing flux; ere the out- 
pide of the core reaches the new value of flux density immedi- 
ately while the center of the core remains at the original value 
and gradually attains the same value specified by the boundary 
in the steady-state condition. <A more complete discussion of 
the eddy-current problem and its relation to the flux distribu- 


tion is given in the next section. 
2.3 Eddy Currents 


Eddy currents are generated by a change of flux in the mag- 
netic material of a core as shown in Fig. 3, and become a very 
important factor in the analysis of the transient behavior of 
high-speed electromagnetic circuits. These currents affect the 
change in flux by tending to produce an opposing flux, thus 
changing the flux distribution across the pole face in addition 
to delaying flux build-up or decay. Flux build-up is not 
affected nearly so significantly as flux decay because the eddy- 


current opposition flux is a small part of the total flux 
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applied. Since a no-work magnet uses flux decay to accomplish 
its function, the eddy currents play a large part in controlling 
the release operation; thus they can limit the speed of opera- 
tion for many of the electromechanical systems used in modern 
input-output machines of digital computation systems. 

One can visualize the effect of eddy currents on the flux 
distribution by considering a cylindrical core to be made up of 
concentric shells, each shell being a hollow cylinder of differ- 
ential thickness. Each shell constitutes a short-circuited turn 
enclosing part of the core flux. The outside shells link all 
the flux while the inner shells link only 6 small part of it. 

A voltage is induced in the shells upon a change of flux. The 
magnitude of the induced emf is determined by the amount of flux 
linked by the coil; thus the outer shells have larger emf's and 
eddy currents than the inner shells. Induced currents apply a 
magnetomotive force to the part of the core that lies within it; 
thus the center of the core is subject to the magnetomotive 
force of all eddy currents while the surface is subject to none. 
This accounts for the distribution present in the transient 
state. Induced currents tend to oppose any change in flux; thus 
they will tend to sustain a decreasing field and oppose an in- 
creasing field. It can be seen that the flux distribution is 
altered by the eddy currents and is not only a function of time 
but also of core radius. An approximate flux distribution with 
radius is shown in Figs. 5 and 6, Figure 5 illustrates the 
effect of eddy currents during the transient flux build-up con- 


dition while Fig. 6 illustrates the effect of eddy currents 


13 B, Pi + Be 


(d) 


Fig. 5.. Effect of eddy currents on flux 
buildup in a ferromagnetic material. 


Fig. 6. Effect of eddy currents on flux 
decay in a ferromagnetic material. 
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during flux decay. 
Let us represent the applied field intensity with vectors 
in the direction of the applied field. We may then illustrate 
the flux at ‘'t = Ot by Fig. 5a. Once the field is applied, 
eddy currents i,j, io, iz, + « «, are set up in the core as 
illustrated by Fig. 5b. These eddy currents induce an opposing 
flux Y, shown in Fig. 5c. The resultant flux at t= 0+ is 
then #) + Z, and is illustrated by Fig. 5d. Since the eddy 
currents depend on a changing flux they will decrease as time 
increases, thus g. will approach zero in the steady state and 
the steady-state flux will be equal to the applied flux g,- 
Upon release the same phenomenon occurs; however, the core 
is originally in the magnetized state and eddy currents will be 
induced opposite in direction to that of the previous example. 
This follows from Lenz's Law. The only flux present at t = 0* 
is that induced by the eddy currents as shown in Fig. 6b. As 
time increases, g. decreases until the flux returns to zero. 
Actually it will return to a residual value g.. 


»Y 


2.4. Flux Distribution in the Rectangular Toroid 


With the preceding discussion of eddy-current phenomena in 
mind, one can readily visuslize a flux distribution existing in 
a toroid of rectangular cross section. Furthermore, one knows 
thet the flux distribution is changing with respect to time and 
values of flux at any point within the core depend on the dis- 


tance from the outer edge of the core at which the boundary 
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condition is applied, the boundary condition itself, the material 
of which the toroid is made, and time. If one considers a rec- 
tangular cross section to be divided into many small squares of 
nearly infinitesimal area, the value of flux density in a given 
Square within the core is a function of time and is different 
for each area. Time dependence of the flux density may be ob- 

+ aia if one can find a relationship sufficient to specify the 
value of magnetic field intensity for this area as a function of 
time. After division into small areas one may then consider 
each area separately as a rectangular toroid with uniform flux 
density provided the area is very small in comparison to the 
total cross-section area. Values of H, B, and % for each area 
can then be calculated for each value of time according to the 
relationships given in section 2.7. One might note that once H 
is determined, it is a simple matter to evaluate B from the 
hysteresis loop. Multiplication of B and the area correspond- 
ing to the value of B calculated will yield a value of flux 
for that particular area. It then becomes possible to calculate 
total values of H, B, and % for the entire cross section. If 
the above. calculations are made for each value of time in the 
Bee ee state, curves relating the time dependence of the flux 
distribution and total flux can be computed for various boundary 
conditions and the effect of parameter changes can be analyzed 
theoretically. Since no literature was found that expressed the 
magnetic field intensity as a function of position and time for 
the case in question, equation (2-28) was derived. It will be 


referred to as the "Hysteretic Diffusion Equation" and its 
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derivation is included in the next section. 


2.5 Derivation of the Hysteretic 
Diffusion Equation 


The derivation presented requires only basic electromagnetic 
theory and simply applies Maxwell's Equations to the special case 
being considered. To refresh the reader's memory, the time var- 
iant Maxwell Equations are listed below in both point and inte- 
gral form. Overbarred symbols designate vector quantities in 
the cartesian co-ordinate system. This system was the most con- 
venient for the rectangular cross section being considered. A 
co-ordinate transformation of the resulting equation could be 


used if different cross-section shapes are studied. 


Maxwell's Equations 


(Point form) (Integral form) 
(2-1) A. VxE = -OB/dt B, gz - di = (aor) [ - aS 
s 
(2-2) A. YVxH =i + OD/dt a. 6H. a =[ T+) - 
s 


OS ee Cee ae 0 


(2-4) A. ve ae BO D- J po 


Consider a toroid with a very large inside diameter in com- 


wil 
rT 
ro 
w 
wi 
Q, 
tal 
TT 


parison to the bar diameter or an infinite bar of ferromagnetic 
material which undergoes a change of magnetomotive force on its 
boundary due to a change of the exciting current I supplied 


from 6 voltage E. (See Fig. 7.) 


Pig eas 


Rectangular toroid. 
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Fig. 7b. Cross section 
of a rectangular toroid. 


Fig. 8. Hysteresis loop. 
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It is known with some certainty that a flux distribution 
exists within the core due to eddy currents while the magnetiza- 
tion of the core is in the transient state, thus stesdy-state 
techniques cannot be applied; however, Maxwell's Equations hold 
in the transient state as well as for the steady state, thus an 
equation relating the effects of H as it varies with x Fy 85 
and t can be derived. 

First we may assume that an mmf has been applied. When 
applied, a magnetic field H, an electric field E, and a flux 
density B, exist in the toroid, and all are governed by Max- 
well's Equations. If the electric field current density D, is 


neglected, Maxwell's second equation becomes 
(2-5) VxH=i 


Since the current flow in a conductor is in the direction 
of the applied electric field and perpendicular to an incremental 
surface as, the total current I is obtained by integrating the 


current density i over the surface; 1.@., 
(2-6) r-fi.B-[@xn.B 
* s s 


These currents are the eddy currents present in the transient 


state and can be represented by the scalar multiplication 
(2-7) i=oEr 


Thus the relationship between the magnetic field intensity H, 
electric field intensity E, material conductivity O , and the 


conduction current density i; is determined since 


18 


(2-8) I 


T 
o— 


i.@- [vx .B-of F.5 
s Ss 


Therefore one can now show that 
(2-9) VxH=OE=i 


A relationship between the flux density B and magnetic 
field intensity H can be found by combining equation (2-5) and 


Maxwell's first equation, equation (2-10). 

(2-10) VxE=- OB/dt 

Combining equation (2-9) with equation (2-10) results in 
(2-11) Vix VY «jn =i-. OB/Ot 


In the general case B, H, and E are functions of Xe 0 ae 
and t and should be written B(x,y,z,t), H(x,y,z,t), and 


E(x,y,z,t). It is also known that a nonlinear relationship 


exists between B and H; thus one can say that 


(2-12) B = f(H) 
cee ae OB/Ot = (OB/OH) (OH/dt) 


Further examination of the functions B(x,y,2Z,t) and 


H(x,y,z,t) illustrates the fact that total derivatives may be 


used since 


(2-1) -— 


dH OH dx OH dy OH dz On 
(2-15) a 
dt Ox at *Oy dt. Os at Ob 
and the derivatives of distance with respect to time are equal 


to zero if the toroid is stationary with respect to the exciting 
mmf. Thus we have 


OB/dt = dB/dt 
OH/Ot = dH/dt 
(2-16) aB/dt 
(2-17) 0B/OH 


(dB/dH) (aH/at) = U(H) dH/dt 


dB/dH = U(H) 


The quantity U(H) represents the relationship between B 
and H as illustrated by the slope of the hysteresis loop of 
Fig. 8. Combining equations (2-17), (2-13), and (2-11), we 


find that the following equation exists. 
(2-18) OH/Ot = -(1/ou(H)) (Vx 9x #) 


Solution of equation (2-18) will provide an expression for 
H as a function of the position within the rectangular toroid 
and time; however, a solution is very hard to obtain. If one 
assumes the only component of field intensity present is along 
the z axis (see Fig. 7) and that this value of field intensity 
H is the same value at all points on the z axis, equation 


(2-18) reduces to a form generally recognized as the diffusion 


equation in two dimensions; i.e., 
(2-19) OH/dt = V* H/o-u(H) 


The above assumptions require that 
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where in general H is given by the vector equation 
Bey Be Hy J tis 


Thus equation (2-18) may be written as 
OH 1 ele te 
ge! eV Ra ee) 


(2200) — 
Ot oO U(H) 


Expanding the expression (Vx V xB) of equation (2-18) with 
H = Hk and noting that 


Vines. Virsa  Gateen be PPG. )H 
Vixn V xe Vil oe ee 
(2-21) Vx VxHe= - Ve Hk 


Thus equation (2-18) is now reduced to the diffusion equation 
for the two-dimensional case, i.e., 
OH 1 O°H, = d°H 


(2-22) at yee + 
Ot oO U(H) Ox? oy? 


An approximate solution of the above equation may be accom- 
plished through use of numerical integration techniques once a 


suitable approximation for U(H) is determined. The remainder of 
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the thesis will deal with the special case outlined sbove. The 
vector notation may be deleted since the only component present 
is H, and the equation was solved only for this particular 
case. | 

The only remaining unknown is the functional relationship 
existing between B and H. One knows this relationship is ex- 
pressed by the hysteresis curves of various magnetic materials; 
thus an equation which approximates the specific curve of inter- 
est would be desirable. Development of such an equation is 


undertaken in the next section. 
2.6 Approximation of the B-H Curve for U(H) 


Because of the flux distribution within the toroidal core, 
esch small area as defined in section 2.4, will have a different 
value of field intensity and flux density at any given time in 
the transient condition; thus an approximation would allow 9g 
to be calculated directly for each area once H for that area 
is known. This may be accomplished if an equation can be found 
which approximates the particular B-H curve of interest. Since 
this particular curve is determined by the initial values of H 
for a given material, an approximation of only this B-H curve 
would be sufficient to compute B once H is known. 

A modified Froelich approximation equation sas given by 
equation (2-23) can be used for generation of one-fourth of a 


given hysteresis loop. 
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C,H 
(2-23) oo 


=——— + B, 
Co + H 


where B is the flux density 

is the residual flux density 

H is the field intensity 

C, constant specifying material properties 


C5 constant specifying material properties. 


Values of C;, and Co determine the shape of the curve 
and may be varied to generate an approximation to most hysteresis 
loops. They may be calculated from B, and selection of two 
values, B,, Hj, and Bo, Ho, taken from near the knee of the de- 
cay portion of a given hysteresis loop representative of the 
operating range in a given material. (See Fig. 8.) 

Since both points selected must lie on the curve, a simple 
simultaneous solution of two equations formed by substituting 
B,, Hj, and By, Ho, in equation (2-23) will be sufficient to 


specify Cy and Co. Performing this operation yields 


C,H 
(o-288) SP Bae ae ag. Bo 
Co tT Hy 
C,H 
(2igaie Wea ee 
CA. Wid 
2 2 


Solving for C, and Co we have the following relationships: 


HoH3(B3 - Bo) 


(2-2h) a 
(By - B,)H3 - (Bz - B,)Ho 
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roceis a (Bo = Bp) (Co + Ho) 
Ha 
After determining C, and Co, one-fourth of the hysteresis 
loop can be plotted using the values calculated from equation 
(2-22). If there is a reasonable correspondence between the 
original curve and the curve plotted using the approximation 
equation, one can assume C, and Co are sufficiently accurate 
to specify the B-H relationship. Reflection and translation of 
this quarter section of the hysteresis loop generates the remasin- 
ing portions of the loop. Calculation of B for values of 
H <-H, for flux decay and H >H, for flux build-up are con- 
sidered in section 4.3. 

This method of approximation is perhaps rather crude; how- 
ever, it does suffice in this case. More accurate approximations 
are no doubt possible although maybe not practical since the mag- 
netic properties of a given material will vary and cause larger 


errors than those due to the approximation. 
2.7 Hysteretic Diffusion Equation 


If we now use the approximation equation (2-22) for 


evaluating U(H), we have 


dB d CH 
(2-26) U(H) Se es cadet + Br 
dH dH |Co +H 


Performing the differentiation we have 


eh. 


CC 
ie 
(2-27) U(H) Aigo. ae wap 
(Cy + H) 
Combining equations (2-22) and (2-27), we obtain the hysteretic 


form of the diffusion equation. 


OH (Cc He) Leen 0°H 
(2-28) z = eae eae 2 ae a 
Ot OC1Co Ox? oy? 


It is this equation whose solution will yield the time 
dependence of the flux distribution pattern in the rectangular 
toroid subject to boundary conditions as outlined in section 2.5. 

The remainder of this thesis is concerned with the numer- 
ical methods for obtaining this solution, implementation of 
these methods, and demonstration of the numerical process by 
actually obtaining a solution of equation (2-28) with its 


boundary values. 
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3.0 NUMERICAL SOLUTIONS WITH THE 
MODIFIED EULER METHOD 


3.1 Euler Method for Solving an 
Ordinary Differential Equation 


Numerous methods of obtaining solutions to ordinary and 
partial differential equations have become practical since the 
advent of the high-speed digital computer. The Modified Euler 
Method was chosen for this problem because of its simplicity, 
although other methods may provide a more accurate approximation 
to the solution or take less computation time. Before showing 
the application of Euler's Method and its modification to 
equation (2-28), let us consider the solution of an ordinary 
differential equation of first order. In symbolic form we may 


write 


(3-1) — = f(t,H) =D 


The integral of equation (3-1) gives H as a function of time; 
thus we have H = F(t). A graph of F(t) is a curve in the 
H-t plane which may be approximated by a series of short line 
segments provided the curve is continuous; thus we have the 
approximation relation (see Fig. 9) 

dH 
(3-2) SH = At tan6@ =|—] At = Dp At 

dt]o 


(3-3) H, = Ho + Do At 


If we let At = h = ty,, - ty, we can express the approximation 
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by 
(3-),) Bee ly OF D;h Bee as Oa ae A re: ge AC) 


This is known as Euler's Method. However, if h is taken small 
enough to yield sufficient accuracy the method is too slow; if 
h is larger, inaccuracies will cause the approximation to be 
unsatisfactory; furthermore, if the graph is monotonic, the 
approximation will diverge from the actual curve for any value 
of h chosen. A modification of this method tends to eliminate 


the divergence. 


3.2 Modified Euler Method for Solving 
an Ordinary Differential Equation 


Starting with an initial value Hj one can approximate 


H, in the same manner as before to yield 
(3-5) H, ‘2) = Hp + Doh 


Substituting H,‘2) into equation (3-1), one obtains an approx- 


imation for dH/dt at the end of the first interval, i.e., 
(3-6). ° p, 2) Cos H,{2)) 


An improved value of H is then found by multiplying h by 
the average of the values of dH/dt at the ends of the interval 


to to ty5 thus we have 


lnote that D; represents (dH/dt),;. This notation will be 
used throughout the remainder of this section. 


27 


Hy 

actual Ho 

actual Hy 
Note: D; = (dH/dt) ; 
ig Pease rs 


t 
to ty i) 
Fig. 9 Approximation with Euler method. 
H 
H, ‘2 eee A ee H = F(t) 
PAL Oe or 
- Paes! AOU Hy 
E 
’ Note: D; = (dH/dt), 
Ho 


Fig. 10. Approximation with modified Euler method. 
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1 


A more accurate value of Hy can now be calculated and will be 


denoted as follows: 


a 


If we look at Fig. 10, it is evident that this value of H, (2) 
is more accurate than H,(), H, {1) is represented by the line 
Ho + BC if calculated according to kuler's formula. Substitution 
of H,{}) in equation (3-1) gives an approximation of the slope 
represented by the tangent at point F. If a value of Hj, were 
calculated using the slope at the end of the interval, we would 
have H) = Ho + BE. When the average of the slopes at the ends 
of the interval are used in place of Do we find that 
Hy‘) = Hp + BG which is definitely a better approximation to 
the real value of Hy, at ty, than the first value of Hj 
computed. This process may be represented symbolically as 

1 


(3-9) AH = — h(Do + Dy!2)) 
}] 


1 vi 
AH = — (BE + BC) = — (BE + BE + EC) 
2 


2 
1 
(3-10) AH = BE + — EC 
(2) Z 
(3-11) H = H, + AH = H, + (BE + — EC) 
iL: O O 5 
lime superscript (k) on H, {«) indicates the k*™® value of Hy 


where H,; is the approximation to the actual value of H, at a 
time t = t; = 0 + At. 
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The new value is much closer to the actual value of 
A, |= (Hp + BF) than before. Continuation of this process by 
again calculating an approximation to the slope at the end point 
of the interval and by substituting H, (2) in equation (3-1), 
will yield a more accurate approximation to the slope at point 


F. A new value of H, 3) can then be calculated. 
(3) “ (2) 


H, ‘3) will be more accurate than H, ‘@) since the approximation 
of the slope will be improved; i.e., p, (2) is more sccurate than 
Dg). Further continuation of the process will yield succes- 
Sively more accurate approximations to the actual value of Hy. 
The process may be continued until the value of ia ~ H,‘*), 
where B, ‘*) is the approximation of the actual value of Hy. It 
must be noted that H, {X) is only an approximation to the actual 
H,. As successive values of H, “) are generated, H,‘*) will 
converge to a value H, {n) which will not be the same as the 
actual value of Hj. To force the approximation H, (*) to con- 
verge to Hy we must make h very small. The approximate solu- 
tion will approach the exact solution as h—®0; however, as 
h-®0O the number of calculations increases and the computation 
time becomes very large. It then becomes necessary to determine 
the magnitude of errors permissible in relation to the time 
available and adjust h accordingly. 

To illustrate the fact that the limit of H, (*) does not 
approach the actual H; as k—® oOo, let us consider the following 


example. Suppose f(t, H) is expressed by equation (3-13) and 
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the initial conditions to = 0.0 and Hp = 1.0. 
(3-13) Ce fh. Mh) 3.) ee 


Substituting these values of tg and Hp in equation (3-13) 


we obtain 


Davart Cigale thei 2 aD 


If we select h 0.05, we may then write 


Pa 32 as 
1) = ol) es 


The second approximation to Hy and Dy are 


= 1.0525 


a8 
bh 
I 
=e 
oO 
t: 
| 
= 
oO 
oO 
4. 
oO 
uh 


oO 
Fe 
ne} 
| 
SF 
a 
aE 
a8 
pa 
I 
ht 
lo 
oO 
ine) 
U1 


Continuing we have successive approximations to Hy and D> 


as follows. 


1 
eds = Hy + - (Dp + p, '2)) = 1.05256 
p,'3) = t, + H,'3) = 1.10256 

1 
ray =H +- h(Do 4 p, '3)) = 1.05256 


2 
Since H,() = H, ‘3? we should ee stop the process if we desire 
H,'*) to agree with H,‘**1) only to the fifth decimal place. 
turther values of H, for 3S k S @Mwill yield the same value 


(x) in the first five decimal positions as H,‘3), This 


(i) 


for Hy 


points out the fact that the approximation Hy will not 


Sas 


converge to the actual H, as ke-®OOsince H, actual = 1.0525). 


We now take 


4, (3) 
p, (3) 


1.0526 
1.1026 


Hy 


Dy 


Continuing, we can calculate the first approximation for Ho 


and Do- 


Ho(2) = H+ D,{2)n= 1.1077 


Do‘1) = t> + Hp ‘2) = 1.2077 


Then we calculate second and third approximations for Hy and 
Dy (i.e., Hp '?) 4, '3) and Dp‘), Dp'3)) and stop since 
Dp‘3) = ba‘*). Consequently Hp 3) agrees with Hp (2) and we 


have 
(3-14) Hy = Ho(3) = 1.1104 
(3-15) Do = Do(3) = 1.2104 


Further continuations will yield successive values of H3, 
D3, Hy, Di,» etc. 

To evaluate the accuracy of the method let us compare the 
approximated values of Hi, Ho» etc., to those calculated from 
the solution to dH/dt = t + H, which is 


(3-16) H @ Deh'=-t <= J 


These are given in Table 1. Accuracies of H indicated in 


Table 1 can be improved only by using a smaller value for h. 
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Table 1. Comparison of approximate and exact solution 
of dH/dt = t+ H=D. 


Oe t : H;act - H; approx : D;act - D; approx 
@) 0.00 1.00000 1.0000 1.00000 1.0000 
Al 0.05 1.0525 1.0526 1.1025) 1.1026 
2 05.10 1.1103) La LO: 1.2103) 1.210) 


In general, execution of the Modified Euler's method for 
the equation D = f(t, H) is as follows. 
Step 1. Obtain the initial conditions 
t = to 
h = constant = At 
H = Ho 
Step 2. Evaluate the approximations to Hj by generating 


the approximations 


(1) (13) (2) (2 (1 k 

D, deE a UdeDy hs HERES Me ) a CEL 
where 
(3-17) Be PD = f(to, Ho) 
(S-1S)h iy Aas) Seg ep 
(3-19) p, **) = (ty, 7 ")) 
(3-20) H, ‘*) =H) +—h (Do + p, ‘*-1)) 

2 

for k= (1, .2,.35 « . «, PD) and p an integer such that 
(3-21) p, ‘P) = p, (Pod) 


(3-22) H, (P+2) = H, ‘P) 


33 


When this condition occurs take 


(3-23) Dy} p, {P-2) 


(3-24) Hy 


( 
Hy Pp) 


Then proceed to Step 3. 
Step 3. Evaluate the approximations to Hj from the 
calculated approximation to H, and D, obtained as the result 


of Step 2 and generation of 


(1) (2) (2) 


Do , Ho ’ Do Pe et) Ho 9 Do 
where 
(k) _ 
(3-25) Hy H, + Dyh 
(3-26) Dok) =. f(t, Hp(X) ) 
z 
(3-27) Hp‘) = Hy + 3 h (Dj) of pa X37) 


for k=1, 2, 3, .. ., q) and q an integer where 


(3-28) Dp'v) a pp‘ 4-2) 
(3-29) Hyatt) = gla) 


When this condition occurs take 


(3-30) Dp = Dp! a2) 


(3-31) Hp = Ho‘ 


Then proceed to Step 4. 
Step 4 will generate the spproximation to H3 and D3 
after which we proceed to Step 5 for an approximation to Hy 


and Dy and so on until we have an approximation for H,, and 


3h. 


th 


D The process should stop when the n approximation has been 


ae 
evaluated. The value chosen for n will depend on the accuracy 
desired and maximum value of At to be used. A flow chart of 
the croceear ae given in Fig. 11 for the function D = f(t, H). 

If the function f(t, H) is specified by the Hysteretic Dif- 
fusion Equation and we apply the Modified Euler Method as before, 


we can solve for the value of H(x, y, z, t) as specified by the 


partial differential equation (2-28). 


3.3 Numerical Solution of the Hysteretic 
Diffusion Equation with the 
Modified Euler Method 


When applying the Modified Euler Method to the Hysteretic 
Diffusion Equation, equation (2-28), we find that evaluating 
successive values of OH/Ot becomes somewhat more complex since 


f(t, H) now depends on many variables; i.e., 


te Se OC1Co Ox? dy? 

Because of this it becomes more difficult to visualize the phys- 
ical significance of OH/Ot in terms of equation variables; 
however, no problem should exist in evaluating OH/Ot if we ex- 
press the equation in its finite difference form and establish 
initial and boundary conditions. The quantity OH/Ot may now be 
evaluated simply by a series of numerical operations. Successive 
values of OH/dt can be found through the recurrence 


t=t-. 
=f 
relations as outlined in the flow chart of the Modified Euler 


START 


ey a - 


~~ 


ayy: 
(ta4a,Hy43 


STOP 


Flow chart of modified Euler method 
for adH/dt = f(t,H) =D. 
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Process. Although some care must be exercised to limit the size 


of At to insure that the 


(3533) a (ot, Se oy eit. 


— 0 
the problem is a straightforward procedure as outlined by the 
fiow chart.of Fig. 11, 
In the following discussion, the magnetic field intensity 
for “the: 7, gth point in the grid of Fig. 12 is specified by 


th approximation to H at the I, gth point is 


th 


Hi dood ):3. thera 


denoted by H(I, J) and the k*® value of the i approximation 
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Ot. Hs Stowe phe kL: gth Pout err) I, D) Phe Note that the super- 
script does not indicate that H(I, J); is raised to a power or 

that 1 is? the ,th derivative of H. Let us also denote (OH/Ot) 5 


by ite 
3.4 Conditions for Convergence 


Conditions exist under which the process will not yield a 


solution. This is due to the fact that lim (H, (et) - H, ‘*)) 
k—&@ 


will not approach zero; however, we can establish a relation- 
ship between the time increment and other parameters such as 
grid size, Ci, Co, Oo jad Hnax to insure convergence of the 
process and existence of a solution. 

Let us divide a rectangular cross section of the toroid 
considered in section 3.0 into a rectangular grid network as 


shown in Fig. 12. Each grid is of equal area h* and H(I, J) 


Fig. 12. 


One-fourth of rectangular toroid cross section for 


sample problem given in section 5.0. (See Fig. 7b.) 


Note: NI =4, NJ = 8, NY = 6, NX = 10. 
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is the field intensity at the intersection of the pe row and 


gJ*2 column of the grid. 


The Hysteretic Diffusion Equation for the T,gth point is 


given by equation (3-3). 


(Co + H(I,5))2 [ O2H(1,g) d2H(I,3) 
ie eae eee i) = ee 


(3-3))+ R(T.) = ss 
O C1Co Oxe oye 


Expressed in its finite difference form the equation becomes 


MCI; De SUB DA) Seated) a HG a) 


UL. oan ud eae, dt+ ea 


Wty) 


(Ooh Mien ees). ¢ V(1, 9) 
———————EE Eee aS 


2 Pe 
O C{Co h h 


(3-35)+ P(I, J) 


If we. wish. to: evaluate 'P(I,J) at l-=°5, J#="3,. we obtain 


V(5,3) = H(6,3) - 2H(5,3) + H(4,3) 
W(5,3) = H(5,4) - 2H(5,3) + H(5,2) 
. (Cp + H(5,3))° [w(5,3)  V(5,3) 
(3-36) Br yay es eee eee ies eee 
oC,C, h2 ne 


3.5 Relationship Between At and Other 
Parameters to Insure Convergence 


To investigate the relationship between At and other param- 


eters of the problem let us consider a sample calculation of the 


lNote that P; represents (OH/Ot),. 
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approximation to H at a corner position of the grid not on 
the boundary, i.e., the position I = NY -1, J=NX-1. (See 
Fig. 7 and Fig. 12, position (I,J) = (5,9).) 

Let us apply initial conditions and boundary conditions 
when all points inside the boundary are at Hp = H, and all points 
on the boundary are at Ho = ~mi,, where m is some real constant. 

If we examine the Hysteretic Diffusion Equation expressed 
in its finite difference form (equation 3-8), we find that the 
largest value of P will occur at the I= NY -1, J = NX -1 
position. To approximate a value of H at this position after 


a time At, we use the equation 


(1) 


(3-37) H(NY - 1, NX- 1), °° = H(NY - 1, NX- 1)p 


+ Po(NY - 1, NX - 1) dt 


Since we are considering only one position in the grid we may 


drop the subscripts; thus 
‘gb ao 
Evaluating Py we find that 


. (Cy Bugz) | O°H - 0°8 


ry +— 
) 
oC,C, Ox? = dy? 
0°H d°x 
since ~ at this position we have 
ox? dy? 
(3-38) Py = 2Xp_ (H(NY , NX - 1) - 2H(NY - 1,NX - 1) 


+ H(NY - 2,NX - 1)) 


ho 


P 


l 


9 = 2Xo(-mH, - 2H, + H,) 
Po — -2X9(m 45 1) Hy 


(3-39) H,'?) = 8, - 2xp(m+1) B, At 


We know the approximation Hy is between -mH, and +H,, thus 


2X_ (mm + 1) Hy, At has to be chosen such that 


We can therefore limit the value of At to satisfy this con- 


G-tions? i Ze, 
(3-41) 0 < 2X, (m tol) Heat Ss (m tol) Be. 


Thus the possible values of At become 


(3-2) OS At S1/2x) 
(Co ‘ a 
where 2 dS a Seca Ty Ar 
O- C1 Coh® 


Let us now investigate a few special cases. 


Case I. ForAt = 1/2x, we have 


Pp, (1) = +2x9(m + 1)H, 


ry 


F4 


ake 


13. Grid section of 


ly. Grid sec 
of H2(9,5)4 
H 


tip” 


of Fig. 12. 
k = 1,2,3, 


i 12. 


Fig. 
of H, adjacent to H,(9,5) 


Indicates values 
at T =, O™. 


(2) 


HW 
aa” 
a 


Indicating values 
i765 DP. 40F 


iH At = 1/2x. 


boundary ~ ~™"p» 


ui 


2 


p,'3) = +2x5(m + 1)H, 


Thus we can see that an oscillation between H, and -mH, 
occurs and the approximation to Hy can never be found since 


the = Jam’ (Hy ee 
k—& ©C 


- Pets pea is not satisfied. 
Case II. Fordt = 1/lx) we have 


Po = -2x9(m + 1)H, . 

H‘2) = Hy - 2x9(m + 1)H, At = = [iran 
p,(1) = 0 

H,'2) = H - xo(m + 1)H, At = f (3 - m) 


ll 

| 
a 

' 
lo 
as 


p,'3) = -x5/2(m + 1)B, : 
H,(4) = Hw -(5/ixo(m + 1)H, At = — (11 - 5m) 


Thus the process is converging to some value Hy §9) and 


lim (H, et) - H,'*))— po is satisfied. One can easily 
k—& CO 


visualize convergence by considering H to be zero on the 
boundary (i.e., m= 0), and establishing approximations H, ‘*? 
as k increases. This process is shown in Fig. 16. 


In order to have convergence we must not include zero or 


1/2Xp in the permissible values; thus 


(3-43) 0 < At <1/2x) 


Selection of the At to be used is dependent on the com- 


putation time available and accuracies desired. As before, large 


+ 


Fig. 


b = Co} 
\c) 
ay : 
H, ‘4 
@) £5 i 
H, (3) 

4. 

“Tip 


15. Grid section of Fig. 12. Indicating values 
of H,(9,5),'*), k = 1,2,3,4, for 


Hooundary = ~“MHp, At = 1/4xp. 


ip 


16. Grid section of Fig. 12. Indicating values 
of Hz (9,5), /, k 1,253, for 


K 
Hpoundary = 0.0, At = 1/4x. 
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At's introduce inaccuracies; however, they may allow evaluation 
of the maximum value of H more rapidly. This is not neces- 
sarily the case since more iterations are necessary to evaluate 
approximations when At is large, and it becomes possible to 
use more time in evaluating these approximations than to evalu- 
ate H in smaller steps of At. Selection of a permissible 
At for convergence requires that 
0-01 Ch? 

(3-4) OSS AD eh ees 

2(Cp5 a Hina x) 

If we force the corner position to converge all other posi- 
tions along the boundary will converge since they only require 
that 

OC, Coh* 


(3-45) OS Beiecr erty came 
(Co . Emax? 


To insure convergence we may write 


G-G, 05h" 


(3-6) 0<=An <= Co Cy Ope ke 
(Co z Bax! 


where the parameters of equation (3-46) are defined as follows: 


Co = constant and is always <0.5 

o~ = material conductivity 

Cy = constant for B-H approximation (see section 
2-6) 

Co = constant for B-H approximation (see section 


2-6) 
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2 
i 


the grid size selected (see Fig. 12) 


Hnax = the magnitude of the maximum value of H(I,J) 
specified by the initial conditions of 
boundary conditions at time t=O. 

It is now evident that At has to be selected in accord- 
ance with the conditions of equation (3-l) and is not an arbi- 
trary variable. It might also be noted that smaller grid sresas 


h2 


will require smaller time increments. This condition 
specifies the accuracy of the method in the same manner 4s 
limiting the value of At to be small in the example of 
section 3.2. 

With the various limitations on At now specified we can 
proceed to implement the solution to equation (2-28) by con- 
structing a program to perform the almost endless number of 


numerical calculations necessary. Discussion of this task is 


given in the next section. 


hé 


4.0 PROGRAM FOR SOLUTION OF THE HYSTERETIC 
DIFFUSION EQUATION 


h.1l Introduction 


The basic function of the program is to execute the Modi- 
fied Euler Method of numerical integration to solve the Hyster- 
etic Diffusion Equation and to obtain the magnetic field in- 
tensity distribution pattern in the cross section of a rectangu- 
lar toroid subject to given initial time and boundary conditions. 
Before this process may be executed we must first define the 
problem in terms of variable names in accordance with the 
FORTRAN IV programming language, read in the input data, and 
establish output data formats. A method of completing this task 
is given in the next section. It is assumed that the reader is 
familiar with the FORTRAN IV language. 

Discussion of the program is given in various sections 
with card numbers indicating locations of the sections in the 
program listing of Appendix C. Notation and symbolism used in 
the program follow closely that used in previous discussion, 
although some modifications were necessary to satisfy program- 


ming requirements. 


4.2 Main Program EMEX 


This program controls all numerical operations and may be 


subdivided to indicate the particular functions performed. 
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.2A Comments (MG90000 - MG90106) 


The cards serve as a partial list of variable names in the 
program and specify various numerical values to provide given 


print-out formats. 


4.2B Input Data (MG90122 - MG90138) 
and (MG90262) 


The first five cards of each data set give all data neces- 
sary to execute a given problem except the initial and boundary 
conditions which are contained on the remaining cards 6 
through n. 

Any variable name beginning with I, J, K, L, M, or N is 
fixed point and read according to an I5 format. Variable names 
beginning with other letters are floating point quantities and 
are read according to an E16.8 format. The E16.8 format sllows 
any of the common systems of units to be used without alters- 
tion of the input-output format specifications. Input variables 


are read in the sequential order as specified below. 


Card 1 EO, HSO, BR, HMAX, CUO 

Card 2 TO, TKM, COND, AX, BY 

Card 3 CDT, KZ, NI, NJ, K, KO, KI, KG 
Card Bes hes. 33) OS 

Card 5 HIO, DHX, I3, Il, I5 

Card 6 

1 ‘ ((HO(I,J), J = 2,NX), I = 2,NY) 
Card n 


line value of n= 5 + (NI + 1)(NJ + 1)/5, NX = NJ + 2, 
NY = NI + 2. 


4.8 


Definition of the variables and methods of obtaining 


numerical values for them are given as follows. 


1. BR 


2 Boy 
BS ys 


3. HIO, DHX 


a. HIO 


4. HMAX 


5. cuo 


6. COND 


4 oie : 


64 BY 


9. NI 


10. NJ 


ll. TO 


Residual flux density specified as the flux den- 
sity for the magnetic field intensity H = 0.0. 


These are values of flux density and magnetic 
field intensity taken near the knee of the decay 
Ses a of the hysteresis curve. (See section 
a a) 


Variables used for control information which 
specify the number of points on the approximate 
hysteresis loop. For a symmetrical hysteresis 
Loop ‘as ‘shown in‘Fig. 22, choose 13 = 

1 + 2(HIO/DHX). 


Specifies the maximum value of H_ to be plotted 
on the hysteresis loop approximation. 


Specifies the increments of H plotted on the 
hysteresis loop approximation. 


Specifies the number of points on the decay por- 
tion of the hysteresis loop approximation. The 
maximum value is 100. 


Specifies the maximum value of field intensity 
applied where H= $/1l=NI/1l. (See equations 
Ee De es BA oes 

Constant to allow for units conversion when com- 
puting values of force. Force = B A/2u0 where 
Ho = 4ncuo. 

Conductivity of ferromagnetic material. 


"X" dimension of the rectangular cross section 
as shown in Fig. 12. 


"Y" dimension of the rectangular cross section 
as shown in Fig. le. 


Number of divisions of width CX in the "X" direc- 
tion. (See Fig. 12, h = CX) 


Number of divisions of width CX in the "Y" direc- 
tion. (See Fig. 12, h = CX.) 


Specifies the initial value of time T(K). 


16. 


17. 


18. 


19. 


20. 


2l. 


. CDT 


KZ 


KI 


KG 


HO(I,J) 
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Subscript indicating iterations. 
as T(K). 


Appears 


Specifies execution for flux decay or flux 
buildup run. (See Appendix C for values.) 


Specifies length of time increment. Same as Co 
of section 3. e. Use CDT < 0.5 to insure con- 
vergence, smaller for higher accurscy. 

Specifies physical position of data set in rela- 
tion to the last set; e.g., if four data sets 
are being run, KZ = 1 for the last set to be run, 
KZ = 4 for the first set to be run, etc. 


Specifies increments of time for which data is 
to be plotted; e.g., if we have 4000 values of 
total force or flux the printer need not plot 
every point; thus it plots every KI point. 
(Use KI such that KO/KI < 800.) 


Specifies increments of time for which a print- 
out of the flux distribution pattern is desired; 
i.e., a print-out for T(0), T(KG), T(2KG), etc. 


Specifies maximum error norm E as given by 
equation (4-1) to be less than EO; i.e., ES EO. 


Specify halt condition. A halt will occur imme- 
ately after any one of the conditions given below 
are met. One should also note that T(K) = K At. 


HS is the sum of H(I,J) for all points in one- 
fourth of the lattice. Halt occurs for HS = HSO, 
buildup; HS S HSO, decay. 


Specifies maximum time T(K) to be considered. 
Halt occurs when T(K) = TKM. 


Specifies maximum number of increments desired. 
Halt occurs when K = KO. 


Specifies data selected for output and format 
specifications to be used. A detailed discussion 
of Ih is given in section 4.2G. 


Specifies numerical value of magnetic field in- 
tensity at every point in one-fourth of the lat- 
tice at T(0) = Ot. This applies the boundary 
conditions and initial conditions at T(0) = 0. 
Quantities are "read" according to the instruc- 
tion ((HO(I,J), J = 2,NX), I = 2,NY). 
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This completes the definition of data set constituents. 
Further information may be obtained by studying the sample 


problem given in section 5.0. 


4.2C Evaluation of the Hysteresis Loop 
Approximation Formula 
(MG901L0 - MG90250) 


This section evaluates the constants C1 and C5 as given 
equations (2-2h) and (2-25), then calculates the B-H relation- 
ship as illustrated by the hysteresis loop of Fig. 22. A 
graphical and printed output of an approximation to the material 
hysteresis loop is available at the beginning of each set of 
output data. To accomplish this task a duplication of the sub- 
routine FLUD was included; however, the variables were changed 
and an output routine was added. Refer to section 4.3 for 


discussion of the procedure. 


4.2D Calculation of Constants 
and Page Heading Routine 


These cards are used to place a heading at the beginning 
of each output data set which defines calculated variables 
specifying operations to be performed; provide a print-out of 
data read from the input data set; check for input data errors; 
etc. Page 76 of section 5.2 was generated by cards contained 


in sections 4.2D and .2E. 
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4.2E Initial Conditions Print-out Routine 
(MG9072u. - MG9087L,) 


This section provides a print-out of the initial input data 
and corresponding initial values of flux density and force at 
T(O) = OF. Cards MG9072) - MG9078 provide the output in a list 
form as shown for Tape 9. The remaining cards give a matrix 
output format as shown for Tape 15. 

In addition to the print-out instructions total values of 
flux (FLUXT) and theoretical force (FORCT) are calculated and 
printed following listing of B(I,J), H(I,J), and FORCE(I,J). 


4.2F Programming the Modified Euler Process 
for the Hysteretic Diffusion Equation 
(MG90876 - MG91060) 


Execution of the process is nearly the same as described 
in section 3.2 and as outlined by the flow chart of Fig. 11; 
however, the condition (P. - (r/r;)) >> 1 (see Fig. 7a) allows 
use of one-fourth of the cross-sectional area because symmet- 
rical boundary conditions are present when the magnetic field 
is applied. 

Tne finite difference representation of the Hysteretic Dif- 
fusion Equation (equation 3-36) requires that H for all points 
in the lattice adjacent to the point under consideration be 
known. We must compute H only st all points inside the outer 
boundary since the boundary conditions require flux on the 
boundary to remain constant. Because of symmetry, values of H 


immediately to the left of the Y-axis are equal to those 
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immediately to the right and values of H immediately below the 
X-axis are equal to those directly above. (See Fig. 12.) Cards 
MG9088) - MG90900 and MG909L.6 - MG90962 were coded to accomplish 
this task when values of I and J placed the point in question on 
the inner boundaries. The desired values of Pp were computed 
for every I and J as were the approximations H, 2), p,{1), 
Hones » e« e, etc. While computing successive approximations to 


H, an accumulative error norm, E was generated. This error is 


expressed by equation (-1). 


NX-1 NY-1 (k+1) (k) 
Wheat ee H(I,J). “BCL J) 
ae sy od | - ‘ 


This method was chosen in order that new approximations for the 
surrounding points would be considered when calculating succes- 
Sive approximations. The maximum error is specified by EO and 
should be small if a high degree of accuracy is desired. Suc- 
cessive approximations for H will continue until the matrix 
error norm, E< KO, When a convergence condition is satisfied, 
NN (variable specifying number of successive approximations to 
H, i.e., the same as (k) in H, *)) is returned to a value of 1; 
values of flux density, total force, and total flux are computed; 
and a new approximation of H = H;, 4 is begun. If the error is 
not converging rapidly enough to yield a solution in a reason- 
able amount of time or if the error is diverging, an error mes- 
sage will alert the operator and a machine halt will occur. A 
flow diagram of the Modified Euler Method for the special prob- 


lem considered is given in Fig. 17. 
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Calculate B-H curve 
Print and plot hysteresis loop 


Calculate constants and initial 
values. Print page heading and 
initial values of variables 


> + Routine A. | 
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Calculate DHO(I,J 
DHO(I, aR = (OH (ar appa 
Po(I,J) = (OnL(T ’3) QT) o 


(DHO(I,J)+DHI(I,J)) 
HO(1I,J)+X7#DT 
E+ABS(T1-H1(I,J) 

Bei so) ae Tl 


HI (IJ) =HO( 1) J)+DHO(L, J) «DT Calculate B1(I 
Calculate Bapees 


] 


\J 


HO (i J) “=H Cais ii) Calculate BL(I,J) ,FORCE(I,J) 
on outer boundaries on boundaries and corners 
) Check for convergence 


Print E, T, NN 


Repeat routine A 


Calculate DELO Led) 
DH1(I,J) = (OH(I,J) OT) 
J) (OT) 


| mee | 
\/ 

<2 

ay 


Gi) Calculate FORCT, FLUXT 


P,(I,J) = (OH(I,J)/OT 
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Print K,FORCT,FLUXT,FORCE(I,J),B1(1I,J),H1(I,J) 


Print values of halt 
variables and denote 
halt condition 


Print and punch values of var- 
iables at last increment of 
time (continuation data) 


on 


Figs 17. 
Flow chart of EMEX. 
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4..2G Output Data (MG91078 - MG9110) 


The remaining cards are related to output data and appro- 
priate headings for identification. Other special control func- 
tions are performed to specify the halt condition and to print- 
out all existing data at halt time. 

Output data are available on Tapes 6, 9, 15, and 16. A 
sample of the output is given for Data Set 3 in section 5.2. 
Each tape has an output data heading for identification and 
prints initial values of H, B, FORCE, total force, total flux, 
and boundary values. Output data contained on the various tapes 
are as given below. 

1. Tape 6 - Contains approximation to hysteresis loop and 
lists total force, total flux, and time. (See 
section 5.2A, Output Tape 6.) 

2. Tape 9 - Lists all values of FORCE, B, and H with the cor- 

responding time. The number of iterations required 
to satisfy the condition for convergence is also 
listed with the associated error norm E, i.e., 
lists values of E, T(K), and NN. 

3. Tape 15 - Provides the same information as Tape 9 except 
that the output appears in a matrix form as shown 
in section 5.2C. Variables to be printed are 
selected by the value of I chosen. (See Table 2.) 

4. Tape 16 - Data on this tape is used to punch a deck which 

may be used as input to a new run continuing from 
the existing values of variables at the completion 
of the first run. A print-out of the same data 
occurs on the output of Tape 6. One must change 
the value of TKM, HSO, and KO to continue. 

Desired output data can be selected by specifying different 


values of Ij. These data are not written on the tapes if not 


desired since the write time is very large when compared to the 
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computation time. I msy range from 1-7 and the output can be 
determined from Table 2. If Ih = 7, only Tape 6 will have a 


useful output. 


Table 2. Values of Ik and related output data. 


Tape 9 (list) : Tape 15 (matrix) 
Ih :H(I,J) :B(I,J) :FORCE(I,J): NN :H(I,J) :B(I,J) :FORCE(I,J) 
1 x x x x x x x 


x x x x 


— Oo WM ew WN 
ral 
ral 


4.3 Subprogram FLUD 


Subprogram FLUD computes a value of B for a corresponding 
H. The approximation is given by equation (2-22); however, the 
value of B is computed by four different methods depending on 
the value of H and whether the flux decay or flux buildup 
problem is being considered. We assume the hysteresis loop to 
be approximated by reflecting the decay portion of the loop and 
translating the entire decay loop to the right by 2H,. (See 
Fig. 18.) Let us consider the following cases to illustrate the 
regions and four different methods of computing B. Figure 18 


illustrates the regions considered. 
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Region 2 B Region 1 


Region 3 Region 
Buildup 


Fig. 18. Approximating the hysteresis loop. 


Fig. 19. Correction of area for 
corners and boundaries. 
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Region 1.7 (Decay run, -H, < H < ©) 


Flux density B is computed by a simple substi- 
tution of H in equation (2-23). 


Region 2, (Decay run, -0O < H < -H,) 
Since Region II of the hysteresis loop is a4 re- 
flection of Region I about the H axis, we may com- 


pute B by evaluating the difference between H 
and H,, forming a new variable equal to the sum of 


-H, and|H - H,| to replace H in equation (2-23). 
Then compute B with H replaced by |H - H,| - H,. 


A minus sign is then assigned to B in order to 
obtain the reflection charscteristic. 


Region 3. (Buildup run, -oO <H< H,) 
Values of B in this region may be obtained by 
adding the quantity 2H, to H and computing the 
negative of B with H replaced by the sum H + 2H. 


Region 4. (Buildup run, H, < H < ©) 
Values of B in this region are computed by sub- 
tracting the quantity 2H, from H and computing 
B with H replaced by the difference H - 2H,. 
Subroutine FLUD determines the particular region in which 
H falls, and computes the corresponding B according to the 


rules discussed above. A flow diagram of the subroutine is 


given in Fig. 20. 
4.4 Subprogram FORCX 


This subroutine computes a fictitious force for each point in 
the lattice where Force = B°A/2i- The reason for calculating 
forces is to obtain a force time relation for the electromechan- 


ical system shown in Fig. 2. The rectangular toroid was used 


lone velue H, is 4@ positive resl number. 
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SUBROUTINE START 


HD. = HC + HD3 


Fig. 20. Flow chart of FLUD. 
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as an approximation to the system with the assumption that the 
gap would not increase the reluctance of the flux path. 
Total force in the toroid is the summstion of forces st 


each grid point; i.e., 
-1 -1 
(4-2) Total force (FORCT) = FORCE(I, J) 
=2 =2 


Computation of the force may be approximated by the formuls given 
if the area A is adjusted as follows. Since approximate 

values of Bare known only for points at intersecting grid 
lines, areas used for computing force values slong all boundaries 
on the lattice should only be one-half of those for the inside 
points and areas used for the corners should be only one-fourth 
of the inside areas. (See Fig. 19.) 

Cards (MG91466 - MG91)80) select the correct formula for 
computing force depending on values of I and J; (i.e., re- 
place A by A/2 on all boundaries and A by A/k for all corners). 

A flow chart of the subroutine is given in Fig. 2l. 


4.5 Subprogram FLUX 


Subroutine FLUX evaluates total flux by forming the product 
B°A. Operations performed closely resemble those of subroutine 
FORCX except only total flux is computed and the flux for esch 
point in the lattice is not. Once B is determined and subrou- 
tine FLUX is called, B is multiplied by the appropriate area 
and added to the existing value of total flux. When all points 


SUBROUTINE START 


FORCE(I,J) = (B1(I,J)°#*DX/2U0) fe 


FORCE(I,J) = (B1(I,J)*#Dx/2U0) /2 @ 


FORCE(I,J) = (B1(I,J)°#Dx/2U0) /y 


RETURN 


END 


Fig. 21. Flow chart of FORCX. 
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have been summed, the total flux is known; however, no matrix 
print-out of the individual flux values can occur since the in- 
dividual values are not stored. A flow diagram would be a dup- 
lication of that for subroutine FORCX with a change of variables. 

A listing of the main program and associated subprograms 
with four sample data sets is given in Appendix C. Further 


information concerning their use will be considered next. 


6h, 
S.0 USING THE PROGRAM ON A SAMPLE PROBLEM 
5.1 Preparing the Input Data 


Perhaps the best method of sR the use of the pro- 
gram is to consider a sample problem of the type outlined in 
section 2.0. Formulation of an input data set specifying all 
input parameters discussed in section 4.2B must be obtained 
before execution. Let us assume we are given a rectangular 
toroid of ferromagnetic material as illustrated by Fig. 7, with 
an initial field intensity of HO. We then reverse the applied 
field such that the boundary at time = Of is at -HO. If HO is 
sufficient to saturate the core a flux reversal will occur and 
a flux distribution will exist in the core during the transient 
state. We can determine the distribution pattern as a function 
of time by solving for H(I,J), B(I,J), and FORCE(I,J) as a func- 
tion of time. Let us assume that we desire to investigate the 
time dependence of the distribution pattern and the relation- 
ship between total flux (FLUXT), total force (FORCT), and time. 
Typical numerical values may be as given below. 

Given: A. Material is 2.5 per cent silicon iron with a hyster- 
esis loop given by Fig. 22. 
2.5 x 107! mhos 


B. Material conductivity (COND) - COND 


C. Material dimensions (see 


Pugs. Le. - AX = 0.00 meter 

BY = 0.0065 meter 

D. Initial conditions - HO = 250 amp-t/m 
TO = 0.0 sec 
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250 amp-t/m 


G. Maximum field intensity - HMAX 
(see section 3.4) . 


H, Boundary conditions 
1. HO(I,NX), I = 
2. HO(NY,J), J = 


1 

= 

oO 
" 


2, NY -250 amp-t/m 
2, NX 

The above information is sufficient to calculate numerical 
values for the remaining input data variables. (See section 
4.2B.) We will use the mks system of units; thus the constant 


bo = kn x 107? webers/amp-t/m = nCU0O; hence 
cuO = 1.0 x 107! weber/amp-t/meter. 


We obtain BR, B2, H2, B3, and H3 from the hysteresis loop 
for the given material. An approximation to this loop is 
shown by Fig. 22. Typical values for 2.5 per cent silicon iron 


are given below. 


BR= 0.71 weber/m* 
H2 = 11.9) amp-t/m 
B2 = 0.80 weber/m* 
H3 = 103.50 amp-t/m 
B3 = 1.20 webers/m@ 


Values of HIO, DHX, and I3 to yield a symmetrical hysteresis 
loop approximation with a range of H from +00 amp-t/m to 
-400 amp-t/m ere 


HIO = 00 
Note: Maximum value of I3 = 100 
DHX = 10 
For symmetrical loop 
13 = 61 I3 = 1 + 2(HIO/DHX) 
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Grid size is selected by specifying NI, the number of divi- 
sions desired in the X-direction. (See Fig. 12.) This also 
specifies the quantity NJ since the grid is square and 
WIS (BY) (AX) /UNT; hence we shall use NI = and NJ = 8. The 
quantity CDT must be less than 0.5 to insure convergence, hence 
CDT = 0.2. Remaining variables of the first five cards control 
the program as discussed in section .2B. We desire to begin 
at T(K) = 0, hence let K = 03 to print-out all data available 
in both formats, hence Ih = 1; to investigate a flux decay, 
hence I5 = 2; to plot every value of FORCT and FLUXT versus time 
T(K), hence KI = 13 to print-out the distribution pattern for 
all variables H(I,J), B(I,J), and FORCE(I,J) for every fifth 
value of T(K), hence KG = 53; and to execute this input data set 
third from the last, hence KZ = 3. Collecting the above 


quantities in a list we have 


CDT = 0.2 K =0 KI =1l 
NI =} Ties, 1 KG = 5 
NI = 8 I5 = 2 KZ = 3 


Variable EO is the maximum error norm E, thus if we wish 
the total error norm for one-fourth of the lattice to be less 
than or equal to 1.0, we require EO < 1.0. Variables HSO, TKM, 
and KO halt execution. (See section 4.2B.) Let us halt when 
T(K) becomes = 20.0 milliseconds, or when K 2100, and not stop 


for the condition HS S HSO, hence we have 
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EO = 1.0 
. Note: 1. T(K) =K At 
HSO = -1.0 x 10 
oO C1 Coh* 
TKM = 20.0 2. At = CDT 3 
, (Cy + HMAX) 
KO = 100 


This completes formulation of a typical input data set. 

If we collect these values, sequence them according to the 
"read" instructions (see section .2B) and express the numbers 
in accordance with the format specifications, we have a data 


set as given by Table 3. 


5.2 Output Data 


Output data are available on Tapes 6, 9, and 15. The out- 
put data set included was generated by execution of the program 
with input data set 3. I = 1 was used to print-out sll saveil- 
able data. Due to the voluminous amount of date available only 
&@ small number of distribution patterns are included; however, 
patterns for all increments of time considered may be obtained 
by specifying KG = 1. Duplications of page headings, etc., were 
cdatads Wilbn the execution is stopped, the halt condition is 
indicated (see page 101) and all data for the last value of 


time are printed-out on all tapes. 


5.2A Output Tape 6 
(Pages 71 - 82) 


Tape 6 lists values for the approximate hysteresis loop 


(pages 71 - 7) and provides a graphical representation of the 


TYPICAL INPUT DATA SET 


DATA SET NO 3 


TABLE 


(5 X 9 MATRIX) 


on PA ER OA OR PA OR 
oo lelelelelelelelole) 
tt ee elie ie ee 
WL WWW OO Ld LW Lu 
(ole) elelelefejeleysls) 
oo elelelelelelelele) 
oo * Talelelelelelefelve, 
oo lelelelelelelelele) 
(ole) COOCO0O0000 
Se leleleleleleleleje) 
oo LOLA LALA LA UV LA USL 
HO NANNNNNANAN | 


ee eeevveveevee 
ee Tees i ee 


ON OO MMAMMMmAMmAN 
C20 G2 BAO0000000 
+1 + $e eeHeeHHH 
wil UW WWW 
CO 0 000030005 
CO 0 CO0000000 
lomo M—ololololololololo) 
C20 0 B0O0000000 
oloml sm olelolololololole) 
olomulocleleiniclolelolo) 
TV NoMo Ml altaltalataltallaltalts 
NS Se nannnaaa 
. 


eeenwee¢*e#e#¢ 
+e ee eee | 


5 


oF SH AMMMAMNMMAMS 
oO ONOOCCOOCCCOCO 
(++ + eee ee ee + 
Wi W WWW 
Soo oo; oooocooce] 
oe .© ESEcececooo 
oO OHQOOOCO0C0000 
SoeeoSo HeSeoeocooeo 
OCOO00 COC0CCOCO000O 
O0OHK0O DBDOOCOCO0CO0O 
KN N NWNNNNNMNnnn 
EN HRNNNNNN NN 
e**O *MOeeeevseeee se 

eee ee +e + | 


Oe NNMMAMMMMMMS 
elol--lelololelelelolelololo) 
+) +e te eteeeeee 
WW Wl 
CO BDDCOODD00000 
olommolololololelelelelole) 
lol solololololololelololo) 
CO CDOD0COCD00000 
olommgolololololololololo) 
olomE lololololelolololole) 
CO SONNNNHKNHAWN 
eae ta latinvlalatia lala in’) 
ee ®@¢,6¢ 6.0 8 4 0.10.89 6 


' Hee eee es | 


HOOOMMMMMAMMMM MS 
ololelelelolelelelelelolelo) 
Seth eee ee eee ee 
WU WW Wd Wb WW Gu WL 
wlelelelelolelolololelelole} 
ainlololelmlolel=lelelelsio) 
Blololololololololelalatelo) 
elololelolelolololelaimiol«) 
ololelolololelelalolelelolo, 
ololelolelelolelololololelo) 
folololololtallaltaltalraltaltaltaltal 
HON DSANAIANANAAAAN 
*eeeseeete#ee#e#*#8e*ee#e¢é¢ 


+H eee eee | 


68 


69 


approximation (page 75). A listing of the total force, total 
flux, and the corresponding time then follows, (pages 78 and 
79). Next is a graphical output illustrating the time depend- 
ence of total force and total flux, (pages 80 and 81, respec- 
tively). At completion of the run, a print-out of the continus- 


tion data available as a punched card output of Tape 16 occurs. 


5.2B pha pila 9 
(Pages 83 - 92) 


The first five pages of print-out for this tape duplicate 
pages 71 - 7h and 76 of Tape 6 output and are not included; 
however, the following page lists initial and boundary condi- 
tions at time T(K) = 0*, (page 83). A listing of corresponding 
values of flux density, force, and time are also given at 
T(K) = O*. As time increases in increments of At, sa complete 
list specifying distribution patterns at T(K) (K = 0, KG, 2KG, 
3KG, . .. ) is listed according to format specifications as 
shown by pages 83 - 92. This tape also contains a listing of 


iterations and their associated error norm E (pages 84 and 86). 


5.2C Output Tape 15 
(Pages 93 - 100) 


Tape 15 gives a more convenient output format since date 
are given in a matrix format which eases interpretation by locat- 
ing values calculated for specific points in the same physical 
location as they would occupy in the matrix of Fig. 12. All 
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data available on Tape 9 except error norm values are available 
on Tape 15. Total values of force and flux are listed follow- 
ing each matrix output. Duplicates of pages 71-7 and 76 of 


Tape 9 precede the output and are not included. 
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5.3 Discussion of Output Data and Results 


For the sample problem considered, a magnetic field in- 
tensity of 250 amp-t/m was sufficient to saturate the core. 
Assuming uniform flux density, initial values of total force and 


flux can be calculated using equations (2-2) and (2-)); hence 


[2es) G = BA = (1.452) (0.00) (0.008) 
= 6.46) x 10-© webers 
Bea (1.4.52) * (0.00) (0.008) 
(Dai) P= —— = = 26. 81,3 newtons 


Total values computed by summing B(I,J)h° and FORCE(I,J) compare 
closely with those given above; thus one may conclude computa- 
tions of total force and flux specified by subroutines FORCX 
and FLUX are valid approximations. 

Distribution patterns are indicated by the print-out of 
Tapes 9 and 15. Tape 15 Output Data indicates values of field 
intensity and flux density for each point in the lattice of 
Fig. 12 in a convenient format. (Note physical position of 
H(I,J) in Fig. 12 and matrix output.) Examining these values 
as time increases, 1.6., for TIME(0O), TIME(S), TIME(10), .. ., 
TIME(100), we note that points near the boundary approach 
boundary conditions much more rapidly than those at the center 
of the core. Given sufficient time, all points will reach the 
same values as that on the outer boundary. Let us think of the 
field intensity at each grid point as a vector pointing in the 


z direction (see Fig. 7), whose magnitude represents the 
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magnitude of field intensity at that point. If s membrane was 
stretched over the tips of the vectors, the surface generated 

as time increases represents the distribution pattern. We may 
visualize this distribution pattern by citing the following 
example. Suppose a soap bubble has been partially formed by 
exerting pressure behind a membrane across a tube of rectangu- 
lar cross section. If we release the pressure sustaining the 
bubble, the membrane will return to its initial state. The 
changing shape of the bubble during the period of time in which 
it is decaying to its initial state, is analogous to the shape 
of the surface representing the flux decay distribution pat- 
terns. We should also note that force depends only on the mag- 
nitude of flux density, thus force is always positive. Further- 
more, it will return to a value corresponding to H = -250 amp-t/m 
when all transients have decayed. 

It is hoped that the solution given for the sample problem 
was useful in illustrating how one might obtain a solution for 
other equations of this type. An explicit expression for H 
and B was not found; however, the thesis does present details 
of how the Modified Euler Method of numerical integration may 
be used to establish a numerical solution which approximates 


the actual solution. 
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APPENDIX A 


Terminology and Formulse 


I. Terminology 


Most symbols are defined in the context of the thesis. 
Definition of variables in the program are listed in section 
4.2B and in the "comments" section of the program listing. 
(Appendix C, Cards MG90000 - MG90106.) Those variables not 
defined or frequently used are given in this section for con- 


venience. Overbarred symbols indicate vector quantities. 


Lis 5 4 k - Unit vectors in cartesian co-ordinates 
207 = = we —) - Del-operator+ 
Ox Oy Oz 

3. E = (Ex, Ey, Ez) - Electric field intensity 

4. B = (B., By, B,) - Flux density 

5. H = (Hy, Hy, Hz) - Magnetic field intensity 
6.D = (D,, Dy» Dz) - Displacement current density 
ao 4 oe (i,, tye i,) - Conduction current density 
8. aS = (dS,,dS,, dS) - Differential surface 

9. dl = (dl,,dl,,d1,) - Differential length 

10. O = - Material conductivity 

ll. Uo = - Permeability of free space 
12. % - Time in seconds 


lPgrentheses denote vector components, e.g., H = (Hy, Hy, Hz) 
= Hi + Hy J + Hk. 


13 
aris 
15. 
Tvs 
cp 
18. 
19. 


203 
eas 
a2 


Zc 


. Curl of 


Bt a) 
H(I,J) 


P,(I,J) 
FORCE(I, J) 
FLUXT 
FORCT 


se 


<q 
~ 
se 
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- Grid size in sections 3.1 and 3.2 

- (dH/dt) 5 

- (3H/dt); 

- Specifies row of matrix shown by Fig. 12 


- Specifies column of matrix shown by 
PigecLes 


- Flux density at grid point I,J 


- Magnetic field intensity at grid 
point I,J 


- OH/Ot at grid point I,J 
- Force at each grid point 
- Total flux 


- Total force. 


II. Formulae 


i j K 
0/ox OPRy dz 
Hy Hy H, 


- Divergence of H 


4k. Other useful relations 


a ¥V-V= 7 

b Vx(V- H) =0 

o. .V «+ ( Vx #) =.0 

d. Vx(¥xH) = (V- #H) - WH 
e. Ax(BxC) - (A xB) xC 
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APPENDIX B 


Units and Conversion Factors 
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Several systems of units used interchangeably are the MKS, 


ccs, 


Length 
Time 


Force 


Voltage 
Current 
Resistance 


Capacitance 


ux 


Flux density 


MMF 


Magnetic 
intensity 


B (lines/cm@) 
B (lines/in*) 
H (oersteds) 


H (amp-t/in) 


and Mixed English. 


MKS 
meter 
second 
newton 
volt 
ampere 
ohm 
farad 


weber 


weber/m* 


amp-turn 


amp-t/m 


Gauss 


Conversion factors are given below. 


x 10° 
2 dw, EEG 
x 103 


x .OOkn 


Conversion Factors 


x 6.4.50 
x 0.155 
x°2.020 


x 0.495 


CGS 
cm 


second 


dyne x 22)8.0x107® 


volt 
ampere 
ohm 


farad 


maxwell 


gauss 


gilbert 


oersted 


MIXED ENGLISH 


CA eS = inch 
x LO = second 
= pound 
> pee 0279) = volt 
page eZ) = ampere 
Kwa 0 = ohm 
> see (Ee = farad 
x 1.0 = maxwell 
x 6.45 = lines/in® 
x 1/.4nx = amp-turn 
x 2.02 = amp-t/in 
Maxwells/in® 


B (lines/in°) 


B (lines/cm?) 


H (amp-t/in) 


H (oersteds) 


Gauss 
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APPENDIX C 


Program Listing with Data Sets 


The program is written in FORTRAN IV language for use with 
the IBM 7090/709 IBSYS operating system. Processing is to be 
with the IBJOB FORTRAN IV compiler. All cards with a "$" sign 
in column 1 are system control cards used for all programs to 
be executed and are not unique to this program. All cards with 
a "C" in columns are comment cards and are not processed. 

The program consists of a main program named EMEX and six 
subprograms FLUX, FORCX, FLUX, PPLT, INSE, and LOGI. The sub- 
programs PPLT, INSE, and LOGI are provided to obtain graphical 
results while FLUD, FORCX, and FLUX are subroutines called by 
EMEX to perform specific calculations. One should also note 
that organization of the program is sequenced by numbers in col- 
umns 76-80 of the listing. These also aided in locating various 
sections within the program in the discussion of section 4.0. 

The data input is on tape unit 5 and data output is on tape 
units 6, 7, 9, 15, and 16. Definition of symbols used is given 
in Appendix A, in the "comments" section of the program listing 
and in section 4.2B. If a larger matrix than 35 x 35 is desired, 
additional storage is required and the dimension statement must 


be altered accordingly. 
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The purpose of this thesis was to investigate the time 
dependence of the magnetic field in a rectangular toroid when 
subjected to an applied magnetomotive force. 

According to Lenz's law, eddy currents generated in the 
transient state induce an mmf in opposition to that applied. 

This results in a nonuniform time varying distribution of field 
intensity and flux within the core during a transition between 
states. An equation describing time dependence of this distribu- 
tion pattern was derived by applying Maxwell's equations with 

the assumption that symmetrical boundary conditions were applied 
along the length of the toroid and no leakage could occur. Com- 
bining this equation with another which provided an approximation 
to the B-H relationship for a given ferromagnetic material, 
yielded the Hysteretic Diffusion Equation. 


OH (Co+H)° [ 02H d°H 


Ot OC1C5 Ox2 Oy? 


Its solution yields the time dependence of the field intensity 
as a function of position in the core. 

The Modified Euler Method of numerical integration was used 
to evaluate time variance of the magnetic field intensity and 
flux distribution patterns from which the time dependence of 
total force and total flux was calculated. 

Discussion of error criteria and selection of an appropri- 
ate time increment to insure convergence were discussed and a 


sample problem illustrating the method was given. 


